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ABSTRACT 

We have used the new 90 GHz MUSTANG camera on the Robert C. Byrd Green Bank Telescope 
(GBT) to map the bright Huygens region of the star-forming region M42 with a resolution of 9" and 
a sensitivity of 2.8mJybeam~ 1 . Ninety GHz is an interesting transition frequency, as MUSTANG 
detects both the free-free emission characteristic of the H II region created by the Trapezium stars, 
normally seen at lower frequencies, and thermal dust emission from the background OMC1 molec- 
ular cloud, normally mapped at higher frequencies. We also present similar data from the 150 GHz 
GISMO camera taken on the IRAM 30 m telescope. This map has 15" resolution. By combining 
the MUSTANG data with 1.4, 8, and 21 GHz radio data from the VLA and GBT, we derive a new 
estimate of the emission measure (EM) averaged electron temperature of T e = 11376 ± 1050 K by an 
original method relating free- free emission intensities at optically thin and optically thick frequencies. 
Combining Infrared Space Ofrservafonz-longwavelength spectrometer (ZSO-LWS) data with our data, 
we derive a new estimate of the dust temperature and spectral emissivity index within the 80" ISO- 
LWS beam toward OrionKL/BN, T d = 42 ± 3K and d = 1-3 ± 0.1. We show that both T d and 
(3d decrease when going from the Hll region and excited OMC1 interface to the denser UV shielded 
part of OMC1 (OrionKL/BN, Orion S). With a model consisting of only free-free and thermal dust 
emission we are able to fit data taken at frequencies from 1.5 GHz to 854 GHz (350 /im). 
Subject headings: Hll regions — ISM:individual(M42, Orion Nebula) — radio continuum:ISM - 
submillimcter 



1. INTRODUCTION 

Th e Orion Nebula (M 42), located 437±19pc from the 
Sun ( |Hirota et al. 2007), is one of the closest regions of 
active high mass star formation. This Hll region, which 
lies in front of the OMC1 molecular cloud, is excited 
by a group of OB stars known as the Trapezium. The 
OMC1 molecular cloud is part of a bigger complex that 
extends over 30° of the sky. It is an ideal site for the 
study of star formation and the physics of the interstel- 
lar medium (ISM). The M42 area has been exte nsively 
mapp ed across the electromagnetic spectrum (see O'Dell 



2001 1 for a review). This area contains hot young stars 
located in the Trapezium, pre-stellar cores, and regions 
of dense molecular gas. With multi- wavelength observa- 
tions the properties of these different regions and how 
they interact can be understood. 

In this paper we report on some of the first scientific 
observations at 90 GHz using the 100 m diameter Robert 
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C. Byrd Green Bank Telescope (GBT). A 9" resolution 
5' x 9' continuum map in the bright Huygens region of 
M42 was made using the new MUSTANG focal plane 
array described in Section [2] To confirm the reliability 
of the MUSTANG images, two independent data anal ysis 
pip eline s were used. These are described in Sections 4.1 
and |4~1 

At 90 GHz the Huygens region is bright, not simply 
with free-free radiation from the ionized gas but also, 
in remarkably equal measure, with thermal dust radia- 
tion from the molecular cloud. However, the two emis- 
sion components have quite different spatial structure. 
We were able to separate them using data at higher and 
lower frequencies between 1.5 GHz and 854 GHz. These 
data are presented in Section [5] We describe the mor- 
phology of OMC1-M42 as seen over this frequency range 
in Section [6] In Section|7j we derive an estimate of the 
electron temperature andin Section[8] we study the dust 
emission. The results are summarized in Section [9l 



2. THE MUSTANG CAMERA 

MUSTAN G, the Multiplexed Squid TES Array at 
Ninety GHz ( picker et al.||2008[ ) is a continuum camera 
built as a user instrument for the GBT. At the heart of 
the instrument is an 8 x 8 focal plane array of Transition 
Edge Sensor (TES) bolometers built at NASA/GSFC. 
Two high density polyethylene lenses reimage the Grego- 
rian focus of the GBT onto the detectors with an effective 
focal length of 162 m so that each detector is spaced on 
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the sky by 4.2", approximately 0.5/A. 10 A 4.2" spacing 
fully samples the sky in a single pointing of the GBT, 
decreasing the need for fast (>0.1Hz) chopping or scan- 
ning in order to reduce 1// noise. Slow movements of 
the telescope (~ 1' per second) produce many redundant 
measurements of each point in the field of view which 
can be used to remove much of the noise on timescales 
from 0.07 s to 0.5 s (the time it takes a source to move 
between pixels and across the array). The scan pattern 
constrains 1// noise on time scales longer than ~5 s (the 
time between observations of the same part of the sky) . 
Lower scanning speeds are a great advantage as large ac- 
celerations can excite oscillations in the GBT's structure 
making accurate pointing problematic. 

MUSTANG has an 81 GHz to 99 GHz bandpass. Its 
optical system provides a uniform illumination of the pri- 
mary mirror of the GBT out to a radius of 45 m and zero 
elsewhere. A best fit Gaussian to the measured beam 
shape has a FWHM of 9", close to the value predicted 
by optical models of the instrument. Below 10 dB the 
measured beam is significantly higher than that expected 
from a perfect telescope. This is due to power being scat- 
tered from small errors in the shape of the GBT's pri- 
mary mirror that occur on length scales of 1-10 m. Dur- 
ing the observations presented in this paper, the surface 
had a 390 /zm RMS, consistent with our measured beam 
efficiency of 10%. Since these observations were made 
there have been significant improvements to the GBT's 
surface and our measured beam efficiency is currently (in 
March 2009) around 20%. 

The output from each detector is obtained at a software 
selectable rate which for these observations was 1kHz. 
The data values are relative to an arbitrary but stable 
zero point and are in arbitrary units of "counts" . The 
conversion of counts to flux units depends on each de- 
tector's gain which varies depending on several factors 
including the bias voltage and the location of the TES 
on its transition. The gain of the detectors can be mea- 
sured using a small black body, "CAL" , located at the 
Lyot stop so it uniformly illuminates all the detectors. 
Tests have shown detector gains to be stable to better 
than a few percent over many hours. 

3. OBSERVATIONS 

The MUSTANG map is the result of 5.6 hours 
of integration time yielding a final map RMS of 
2.8 mJy beam -1 . Observations were conducted over 4 
sessions in January and February 2008. Each observing 
session was begun by mapping the primar y calibrators 
Mars or Saturn. In this paper we use the |Ulich| ( |1981[ ) 
measurement for the 90 GHz brightness temperature of 
Saturn, Tb = 149.3 K. Every 30 minutes focus and point- 
ing corrections were measured and applied by collecting 
small maps of a bright source at a range of focus settings. 
In cases where these data indicated that the pointing cor- 
rection was significantly less than MUSTANG'S instanta- 
neous field of view, corrections were applied in the data 
analysis offline. The bright quasar 0530+135 was used 
for this purpose as it is unresolved and is located only sev- 
eral degrees away from M42. The gains of the detectors 
were measured by taking a 30 second scan while pulsing 

10 During the summer of 2008 this spacing was increased to 5.7" 
in order to obtain better signal to noise. 
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Fig. 1. — The GBT boresight trajectory on the sky during the 
"box" scan pattern used in our observations. The square in the 
bottom left represents MUSTANG'S 33.6" square field of view. 

CAL at 0.5 Hz. Another 30 second scan on a blank piece 
of sky was taken to aid with estimations of the weights 
used when co-adding data from different detectors. 

When scanned across the sky a fully-sampled imaging 
detector array such as MUSTANG produces many re- 
dundant measurements of a given Nyquist sky pixel. In 
addition to the sky signal, systematic signals (e.g., at- 
mospheric emission and gain drifts) are present. For a 
well-designed system these will change slowly compared 
to the sky signal and have a different characteristic sig- 
nature in the detector time stream. This is facilitated 
by an appropriate sky modulation (scan) pattern. Fur- 
thermore, most of the dominant systematic signals — in 
particular atmospheric emission and thermal drifts in the 
receiver — are common-mode between detectors. Pro- 
vided that the sky map is highly oversampled on most 
spatial scales it is, in principle, straightforward to form 
high-quality astronomical images. 

M42 was mapped primarily with tiled 5' x 5' "box" 
scans which take 5 minutes each. The telescope is 
scanned diagonally across the rectangular area being 
mapped at a constant speed, turning around at the edges 
with a gentle acceleration so as to minimize resonances in 
the GBT's structure (Figure 111. The scan speed is cho- 
sen so that all astronomical objects with spatial scales 
smaller than a few arcminutes pass through the beam 
in under 5 seconds, a timescale faster than changes in 
the atmosphere or instrumental response. This scan pat- 
tern generates uniform coverage with a high degree of 
cross linking on many different timescales enabling low 
frequency instrumental and atmospheric signals to be re- 
moved. The 5' x 9' selected region of M42 was mapped 
in sets of 6 box scans on subdivided north and south 
fields. To provide additional cross linking, the central 
pointing of each scan was offset from the previous one by 
a non-integral multiple of the spacing between crossings. 

4. MUSTANG DATA REDUCTION 
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Fig. 2. — The M42 MUSTANG map reduced with Obit (Section |4.2| , The large dynamic range, over a factor of 1000, can be seen by 
comparing the linear scale (left) and the logarithmic scale (right). Shown on each map are the Trapezium stars with the bigger one being 
O 1 COri, the main ionization source. 



MUSTANG'S data were reduced using two indepen- 
dent pipelines. The first, written in LDL , implements a 
least-squares map-making approach explicitly solving for 
instrumental paramet ers and sky p ixel values. The sec- 
ond, known as Obit 12 ( Cotton|2008 1, is an iterative algo- 
rithm that uses a single-dish implementation of CLEAN. 
The following subsections describe each approach. Sev- 
eral general properties of the MUSTANG data are worth 
noting: 

1 . The largest signal component is due to atmospheric 
emission. Due to the large size of the GBT and the 
small divergence of the beams, each detector sam- 
ples the same volume of atmosphere. Once the de- 
tectors are properly calibrated this signal is almost 
completely common mode. 

2. There is also a strong component due to MUS- 
TANG'S pulse tube cooler causing temperature 
fluctuations in the optical path. This signal is 
nearly perfectly sinusoidal and its frequency is very 
stable at 1.4117±0.005Hz. The best way to remove 
this signal was found to be by fitting a sine wave 
with a slowly varying amplitude to the low pass 
filtered time streams. 

3. There are significant (0.1 — 0.2 second) absolute 
timing offsets in the MUSTANG data. At our typ- 
ical slew speeds of 1' per second this corresponds 
to pointing offsets in the scan direction compara- 
ble to the MUSTANG beam. A 1 pulse-per-second 
(1PPS) signal is available in the GBT receiver room 

11 http://www.ittvis.com/ProductServices/IDL.aspx 

12 http://www.cv.nrao.edu/ bcotton/Obit.html 



which demarks UTC second boundaries. This is 
sampled synchronously with the detectors and used 
to correct the instrument's timestamps offline. 

Since most of the systematic effects are common mode 
and a naive common mode subtraction will remove any 
sky structures larger than the instantaneous field of view, 
we wanted to carefully test the efficiency of our imaging 
algorithms at retrieving the large scale emission. Conse- 
quently we analyzed the data with two independent ap- 
proaches and compared the resulting maps, as described 
in the following sections. 

4.1. Least Squares Map Making 

The map-making problem can be framed as a large 
linear least squares problem. As an example consider 
that we have Ndet detectors, each yielding N int integra- 
tions for a total of N to t = Ndet x Ni nt data points. The 
time stream data can be represented by a single data 
vector A with N to t components. The observation cov- 
ers a map comprising N P i X <C N tot sky pixels. Further- 
more we allow N par N tot parameters to describe in- 
strumental systematics — for instance, a slowly varying 
common mode due to sky emission and a set of fixed 
detector gains. The sky plus instrument model param- 



eters are represented by a vector, S, with N, 



pix 



N, 



par 



parameters. If the instrument model is linear, or lin- 
earizable, the data are related to the model parameters 
by a Ntot x (N p i X + N par ) matrix H by the equation 



As long as N tot > N pix 
solution of the form 



A = HS (1) 
N par , and H is non-singular, a 



S = (HWH)- 1 HW A. 



(2) 
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exists where we have introduced a weight matrix, W, 
describing the noise and correlations in the time stream 
data. An advantage of this approach is that it is possible 
to explicitly determine the uncertainties in the resulting 
instrument and sky model parameters as well as their 
correlations. This approach has been widely used for 
imaging the Microwave Backgro und and in other deep 



extragalactic exp eriments (e.g., Tegmark 1997 Patan- 
chon et al.|2008 1 . Two difficulties are present. First, the 
matrix HWH is of size (N p i X + N par ) 2 , which can be un- 
wieldy if treated with brute force. Second, if we desire to 
solve for detector gains, the problem is no longer linear 
in the model parameters, but rather is bilinear due to 
g ain x sky pixel term s. 

Fixsen et _al.| ( |2000[ ) have implemented an iterative solu- 
tion to Eq. |5]which, for typical data sets, can be executed 
quickly on modern desktop computers. We have inte- 
grated their least squares imaging routines into an IDL 
data analysis pipeline initially developed by the project 
team for quick-look and explor atory data reduction . Sev- 
eral models are provided by the Fixsen et al.| ( |2000| code. 
We have explored variants of a model of the form 



A = G(S + A) + F, 



(3) 



where G is a set of detector gains, S is the sky (as- 
tronomy) model, A is a common mode optical loading 
term, and F is a per-detector additive offset. Care must 
be taken that the number of model parameters do not 
become excessive, which can be done by a combina- 
tion of selecting a long solution interval for the variable 
terms, and by fixing or eliminating degenerate terms. For 
the MUSTANG data presented in this paper a compre- 
hensive noise characterization was not necessary. The 
weight matrix W was assumed to be diagonal with per- 
detector variances determined every half hour from the 
time stream model residuals. Imaging fainter extended 
sources will benefit from a more accurate treatment. 

The final least squares map was produced by imag- 
ing roughly half-hour long sets of scans and coadding 
the resulting images, properly accounting for the pixel 
weight distributions. The G x A models have an exact 
degeneracy which corresponds to an arbitrarily oriented 
two-dimensional linear ramp in pixel intensities across 
the map. This was eliminated by fitting a plane to two 
orthogonal, approximately emission- free lines of pixels in 
each half-hour map and subtracting it before coaddition. 

4.2. The iterative CLEAN algorithm — Obit 

The basic approach used in the Obit reduction of MUS- 
TANG data is to iteratively estimate the background sig- 
nal and the astronomical signal (the "sky model"). At 
the start of each iteration loop, the current sky model 
is subtracted from the time stream data and the back- 
ground signals (atmosphere, detector, etc.) are esti- 
mated by low-pass filtering the residual. These back- 
grounds are then removed from the original time streams, 
an observed sky image is produced, and the true sky es- 
timated using CLEAN. The sky model and the back- 
grounds are iteratively refined by reducing the timescale 
of the low-pass filter used in the background estimation 
as the iteration progresses. 

Time stream data are imaged by multiplying each da- 
tum with a griding kernel producing a continuous func- 
tion which is then sampled onto a regular celestial grid. 



Two grids are accumulated. The first grid, the image, 
contains the data samples multiplied by a data weight 
and the griding function. The second grid contains the 
data weights multiplied by the griding function. When 
all data are accumulated, the image is normalized by di- 
viding by the second grid on a pixel- by-pixel basis. The 
griding kernel was chosen to minimize the broadening of 
the instrumental response and is an exponential times a 
sine function. 

There are two types of operations used by Obit, those 
done independently on each 5 minute scan and those 
which are coupled to all data through the sky model. 
The independent processing steps are the following. 

1. Reading archive data. The kHz sampled data are 
read and averaged to a 20 Hz rate. Then the 
1.4117 Hz signal from the pulse tube refrigerator is 
estimated and removed by fitting a sine wave with 
a common phase and a detector-dependent ampli- 
tude. 

2. Gain calibration. Detector gains are determined 
from the half hourly CAL scans. These gains are 
applied to all subsequent data until the next cal- 
ibration scan. Calibration to Jansky is achieved 
from observations of planets. 

3. Weight calibration. Statistical weights are deter- 
mined for each detector during blank sky observa- 
tions by the RMS residuals following the fitting and 
removal of a low order polynomial from the data. 
Dead and badly behaving detectors are given zero 
weight. 

4. Common atmosphere and detector offsets. A single 
detector offset per scan and a time variable com- 
mon "atmospheric" term is estimated by low pass 
filtering the median of all working detectors in each 
time sample. 

When these operations are complete, multiple scans 
over multiple days can be jointly imaged and non- 
astronomical backgrounds removed. The imag- 
ing/background removal cycle consists of the following. 

1. Imaging. The current estimates of the backgrounds 
are removed and all data are accumulated onto a 
common grid and normalized to form the "Raw" 
image. 

2. Deconvolution. CLEAN is used to generate a sky 
model consisting of delta functions. Fixed or inter- 
active specification of the CLEAN windows is used 
to guide the process. 

3. Subtracting the model. The CLEAN delta func- 
tions are convolved with a Gaussian approxima- 
tion to the telescope beam shape. The value corre- 
sponding to each data sample is interpolated from 
this CLEAN image and subtracted from the time 
stream data. This results in a residual time stream 
in which the current estimate of the astronomical 
signal has been removed. 
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4. Estimation of the background. An estimation of 
the background signal is made by low pass filter- 
ing of the residual time stream data. This filtering 
alternates between estimating individual detector 
offsets and a common mode offset estimated with 
shortest time scales one third of that used for the 
detector offsets. 

After the cycles have converged, the final CLEAN im- 
age is used for subsequent analysis. 

4.3. Comparison of Results 

The map produced by the Obit pipeline is shown in 
Figure [2] On a linear scale the Obit maps are very 
similar to those produced by least squares. The peak 
height of Orion KL, the brightest feature in the maps 
was within 5% of each other and neither data pipeline 
produced significant negative flux. For bright features, 
all extended structure on angular scales up to 2' was 
recovered. The most prominent discrepancies between 
the maps were faint extended structures in the north of 
the maps which have higher signal-to-noise in the Obit 
maps. Given the similarity of the maps produced by in- 
dependent analysis we believe all of the features above 
the stated noise level to be real. 

5. ORION AT OTHER FREQUENCIES 

To interpret the MUSTANG data we made compar- 
isons to maps at other frequencies. These included 
a 854 GHz (350 H SHARC-1I map and the 664 GHz 
(450 Mm) and 352 GHz ( 850 Mm) maps from SCUBA 
( jjohnstone fc Bally|[l999| . W e obtained new 150 GHz 
maps using data from GISMO ( |Staguhn et al.|2008j and 
radio frequency data at 21.5 GHz, 8.4 GHz, and 1.5 GHz 
(K band, X band, and L band) from the GBT and the 
VLA. A summary of the angular resolution, RMS, and 
calibration uncertainties of each of these maps is given 
in Table [T] In order to compare them, all the maps were 
convertea to MJysr -1 by dividing by the solid angles of 
their respective beams. 

5.1. The 1.5 GHz map 

The image of M42 at 1.5 GHz (20 cm) was derived from 
archival VLA data taken with the VLA in "B", "C", 
and "D" configurations. This combination of configu- 
rations gives sensitivity to angular scales between 14' 
and 4". Initial calibration of the data were carried out 
in AIPS 13 using standard VLA flux density calibrators 
(3C48, 3C138, and 3C2860) and B0500+019 as an astro- 
metric calibrator. 

After the individual observations were calibrated, they 
were concatenated and processed in Obit. Automated 
flagging of RFI was followed by a n imaging using the 



"SDI" technique ( |Steer et al.| 1984 1 . This technique uses 
a standard visibility- based CLEAN to the point that the 
fine scale structure has been removed from the residual 
and only a large, flat plateau is left. In subsequent iter- 
ations, a CLEAN component is added at the location of 
each residual pixel above a given limit. This suppresses 
instabilities inherent with CLEAN that can produce arti- 
facts for very extended sources. Twenty one facets were 

13 http://www.aips.nrao.edu 



used to cover the full field of view as well as outlying 
sources. Twenty million CLEAN components with a to- 
tal flux density of 324 Jy were used in the deconvolution. 
A Briggs Robust weighting parameter of -0.5 gave an ap- 
proximately 7" FWHM beam and the restoration used a 
7" FWHM Gaussian beam. 

5.2. The 84 GHz Map 

The 8.4 GHz image was made by combining continuum 
observations from the GBT and the VLA. Data from the 
GBT was taken in an "On-The-Fly" raster mode and 
contains information on angular scales from the size of 
the map down to the limit of the GBT's resolution at 
8.4 GHz (1.46'). The VLA data, taken with the array in 
"D" configuration, are a series of tiled 3' x 3' pointings 
and contain information on spatial scales from those of 
the shortest to the longest baselines (2' to 8.4"). These 
data sets were combined using AIPS++ and CASA 14 to 
produce a single map with a very high dynamic range and 
sensitivity to all structures down to a FWHM resolution 
limit of 8.4". 

5.3. The 21.5 GHz Map 

A 21.5 GHz continuum map of M42 was obtained from 
archival GBT data and reduced in IDL. The data were 
collected using the GBT's Digital Continuum Receiver 
(DCR), a general- purpose continuum back end. An in- 
ternal noise calibration source with a known noise tem- 
perature was pulsed once within each integration and 
the data were calibrated to antenna temperature using 
this noise source. Data were collected with simple, cross- 
linked raster scans 12' long in Right Ascension and Dec- 
lination. A DC offset and a gradient were removed from 
each scan, and the resulting data gridded onto the sky. 
The resulting map contains information on angular scales 
from 6' to 33.5", the resolution of the GBT at 21.5 GHz. 

5.4. The 150 GHz GISMO Map 

GISMO is another TES bolometer camera bu ilt to op- 
erate at 150 GHz on the 30 m IRAM telescope ( Staguhn 
et al. 20081. The data presented in this paper were ob- 
tained during the first run with the instrument in Novem- 
ber 2007. A "data cleaning" routine was used to remove 
spikes in the data and to identify bad pixels. Since the 
weather was exceptionally good (r 2 25 = 0.07) a high-pass 
filter with a low frequency cutoff at 0.2 Hz removed most 
of the atmospheric noise. After this step a visual inspec- 
tion of the data of each individual pixel was carried out 
and pixels with spikes or steps that were not already iden- 
tified by software were flagged as bad. Per-pixel gains 
and offsets were then calculated by fitting a linear func- 
tion to the correlation between the time series of data for 
each single pixel and the time series of median values for 
each frame. Pointing and flux calibrations were obtained 
by observations of quasars and planets. The total inte- 
gration time for the part of the map shown in this paper 
was about 4 minutes. A combination of fast scanning 
and GISMO 's 2' x 4' instantaneous field of view ensured 
that all extended structure up to (and probably beyond) 
an angular size of 2' was recovered. 

14 Information on CASA - The Common Astron- 
omy Software Applications package is available at 
http: / / science.nrao.edu / index. php/Data %20Processing/CASA 
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TABLE 1 

A SUMMARY OF THE DIFFERENT FREQUENCY MAPS USED IN THIS PAPER. 



Map 


A 


V 


FWHM 


largest Angular 


Relative 


RMS 


RMS 




(mm) 


(GHz) 


(") 


Scale (') 


uncertainty 


(MJysr- 1 ) 


(mjy beam - J ) 


SHARC-II 350 


0.35 


854 


9 


> 2 


30% 


509 


1100 


SCUBA 450 


0.45 


664 


8 a 


1.1 


30% 


106 


300 


SCUBA 850 


0.85 


352 


15° 


1.1 


15% 


5.8 


40 


GISMO 


2.0 


150 


15 


> 2 


20% 


17 


100 


MUSTANG 


3.3 


90 


9 


2 


15% 


1.3 


2.S 


K band 


11 


21.5 


33.5 


6 


15% 


6 


180 


X band 


35 


8.4 


8.4 


2 


15% 


3.7 


7 


L band 


200 


1.5 


7 


14 


15% 


2.8 


3.6 


3A beam shapes consist of a core Gaussian with a FWHM given in this table and a broad component with a 
intensity of 5% of the core and a FWHM of 30" l|Johnstone et al.|2006}. 



5.5. The SCUBA Maps 

We used maps of the Orio n ridge originally observed 
by Johnstone fc Bally] ( |1999[ ). The data were taken of a 
50 long region over four nights with the SCUBA instru- 
ment on JCMT. The observations consist of a number 
of 10' by 10' maps made with chop throws of 20", 30", 
and 65". By combining maps made with these different 
chop throws and different chop directions all extended 
structures with angular sizes up to the largest chop size 
(65") were recovered. 

5.6. The SHARC-II Map 



We used the SHARC-II camera ( powell et al.||2003 l 
on the Caltech Submillimcter Observatory to obtain a 
350/im (854 GHz) map. The absolute uncertainty for this 
measurement is ~ 30% and the beam size is 9". SHARC- 
II's field of view is 2.59' x 0.97' so it can be expected that 
all extended structure with angular scales less than 2' is 
included in our map. 

6. THE MORPHOLOGY OF M42-OMC1 

Figure[3] displays maps of the eight datasets we used. 
Each dataset is presented at its nominal resolution. At 
21.5, 8.4 and 1.5 GHz (K band, X band, L band), the ob- 
served intensity is entirely dominated by free-free emis- 
sion from the ionized gas around the OB stars of the 
Trapezium, since there is no evidence for significant syn- 
chrotr o n emi ssion in M42 as noted by Subrahmanyan 
et al.| ( |2001[ ). At higher frequencies, the SHARC-II 
r^5TGHz7and SCUBA (664 and 352 GHz) maps are 
dominated by thermal dust emission of the OMC1 molec- 
ular cloud lying behind the ionized gas. The two bright- 
est features seen with the thermal dust emission are 
Orion KL/BN and Orion S. At 90 GHz, free-free and 
thermal dust emission have about the same intensity in 
this region, and so both the ionized gas and the OMC1 
molecular cloud are seen in the MUSTANG map. As the 
spectral shape of the dust emission is quite steep com- 
pared to free-free emission, this equipartition vanishes 
rapidly when looking at a different frequency. Thus, the 
spatial structure in the GISMO map is dominated by the 
OMC1 features with only a relatively low contribution 
from ionized gas. 

To the southeast of the Trapezium stars is the well 
known "Orion bar" associated with the edge-on bound- 
ary between the Hll region and the molecular cloud. A 
close up of the bar is shown in Figure [4] The clear shift 
in the position of the bar at different frequencies is what 
one would expect at the ionization front. Moving away 



TABLE 2 

a comparison between our em valu es and those of 
[Baldwin et al.| ( [1991] ). 



location relative 


Baldwin et al. 


this paper 


to 9 1 On C 


(pc cm -6 ) 


(pc cm -6 ) 


30" west 


8.8 X 10 6 


8.2 ± 1.1 X 10 6 


50" west 


4.3 x 10 6 


3.8 ±0.5 x 10 6 


100" west 


8.8 x 10 5 


5.3 ±0.7 x 10 5 



from the exciting stars, one sees first the ionization front 
emitting strongly at 90 GHz due to high emission mea- 
sure (EM), next comes the mid-infrared emission bands 
of polycyclic aromatic hydrocarbons fro m the photodis- 
socia tion region seen with IRAC 8 /im (Megeath et al. 
2005 1 and finally, the molecular material which is bright 
aTthe SCUBA 352 GHz (850 /im) frequency due to the 
high column density of dust. Note that another much 
fainter bar can be seen in the MUSTANG map to the 
north-east of the Trapezium at 02000 ~ 5 h 35' 23" and 

£2000 5° 20' 20" (Figure^. This bar has already been 

detected at radio frequ encies | Yuscf-Zadch 1990 J and 
with molecular tracers ( Fuente~et al.||1996| |Rodriguez- 
Franco et aL|l998[) . 

T he spectra in Figure|5] are those observed at the free- 
free emission peak, and at the Orion KL/BN and Orion S 
locations. These are the averaged spectra over the ar- 
eas delineated by circles on Figure|3] Before averaging 
each map was brought to the spatial resolution of the 
21.5 GHz map (33.5") by convolving it with a Gaussian 
of appropriate FWHM. Spectra were fitted using the 
function described in the Appendix. As noted previously, 
at the MUSTANG frequency (90 GHz), the emission to- 
ward Orion KL/BN has almost equal contributions from 
thermal dust and free-free emission. However, toward 
the free-free emission peak, free-free emission dominates 
the MUSTANG intensity and even the GISMO one at 
150 GHz. 

Toward the free-free emission peak, we obtained an 



emission measure (EM) value of 7.6 ± 0.8 x 10 6 peem 



in a 33.5" FWHM beam 
with a value of 2.7 x 10 6 



pc cm 



10.55 GHz by [Subrahmanyan et al. ( |2001l ). 
account the beam dilution from 33.5 to 90". 
a value of 5.1 ± 0.5 x 10 6 pccm -6 , a factor 
the 



This value can be compared 
6 in a 90" FW HM beam at 
Taking into 
we obtained 
~ 1.9 above 



Subrahmanyan et al. 



(20011) value. This discrepancy 
. d, in part, come from the fact that they used a dif- 
ferencing receiver which could potentially miss flux on 
spatial scales similar to the separation of the feeds. We 



can also compare our EM estimate with those of Bald 



SCUBA 850 



I 

352 GHz 
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Fig. 3. — Maps of M42. All datascts arc displayed with their nominal resolution. The inset boxes in each map display Gaussian 



approximations to the point-spread-functions. The 30" circles on each map delineate areas used to compute average spectra at Orion KL 
Orion S and the free-free emission peak locations (shown in Figure^. Shown on the MUSTANG map arc the Trapezium stars with the 
bigger one being O 1 C Ori, the main ionization source 



mea- 



et al. (19911 that were deduced from HI 1-3 
surements. We derived an EM map (Figure |6| with 8.4' 
resolution using only the 8.4 GHz (X band) and 1.5 GHz 
(L band) data which can be directly compared wi th the 



East-West cut tow ard the West of 9 1 Ori C given by Bald- 
win "eTaTl ( fl99l| ). As can be seen in Table [2] botffThe 
absolute value and the spatial variation of our radio es- 
timate of EM arc in good agreement with their optical 



one. 

We found that no anomalous emission component is 
needed to reproduce data at all frequencies and angular 
scales (< 1.1') considered in this paper. A number of de- 
tections of such emis sion have been made in other parts of 
the sky, for example Kogut et al. ( |1996[ ) found an excess 
between 30 GHz and 90 GHz at higETatitudes in COBE 
data. This excess is correlated with DIRBE far-infrared 
data and has since been interpreted as arising from spin- 
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Fig. 4. — A false color map of the Orion bar. The SCUBA 850 fim image is shown in blue, the IRAC 8 fim map is shown in green, and 
red is the MUSTANG data. The profiles displayed correspond to the cut shown on the map, starting from the upper right and moving 
away from the Trapezium stars which are located toward the northwest of this map. 
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Fig. 5. — Average spectra over the areas delineated by the 30" 
diameter circles on Figure[3] The dot-dashed line is the free-free 
emission component, the dashed line is the thermal dust emission 
component and solid line is the total. For comparison, the dotted 
line shows the free-free emission without taking account of the 
opacity. This fitted function is described in the Appendix. 



ning dust (e.g., Lazarian & Finkbeiner 2003). Excess 
emission on large spatial scales has also been detected 
at fre quencies of 15 GHz and 30 GH z by Leitch et al. 
(|1997|) and at 33 GHz and 43 GHz by |de Oliveira-Costa 
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Fig. 6. — A map of the emission measure derived from the 
8.4 GHz and 1.5 GHz data at 8.4" resolution. The stars show the 
location of the Trapezium stars. 



et al. (19971. Anomalous emission was also observ ed to- 



inson et al. 



200 



2006 2009), cold dense clouds (e.g 



sus et a£]2U06p and planetary nebula (e.g., 





Gas as- 


Casassus 



et al 

peak between 10 GHz and 50 GH z so a small component 
could be missed in our analysis. (Lazarian & Finkbeiner 
2003 1 . In order to check for the presence of anomalous 



emission, new data at Ka band (27 GHz-42 GHz) have 
been obtained and analysis is currently underway. 



7. 



T E DETERMINATION FROM THE FREE-FREE 
INTENSITY RELATIONSHIP 



Assuming that free-free emission is optically thin at 
frequency v\ and potentially optically thick at frequency 
V2, the relationship between intensities at these two fre- 
quencies is 



Iu2 



gff(T e ,is 2 ) 
9ff{T e ,vi) 



h"2 
k T, 



1 



— e 



r(I Vl ,T e< v 2 ) 

,T e ,u 2 ) 



(4) 



wit h Qff (T e , i/ 2 ) being the Gaunt factor defined by Equa- 
tion A5 and j{I y, i T e , U2) the optical depth as defined 
by Equation A6 (with EM oc I Vl ). The top panel of 



Figure [7] shows this function for 1.5 GHz (I V2 ) versus 
21.5 GHz (I Vl ) at different electron temperatures. Fit- 
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Fig. 7. — Modeled (top panel) and observed (bottom panel) re- 
lationships between 1.5 GHz and 21.5 GHz free- free emission inten- 
sities. On the bottom panel, in gray, is the relationship for the 
fitted electron temperature of T e = 11376 K. For comparison, the 
curved dashed line is the expected relationship for T e = 7865 K 
( Subrahmanyan 1992 ) and the straight dashed line is the expected 
relationship without opacity (i.e. T e = oo). 



ting this function to the observed relationship provides 
T e , formally the only independent parameter of Equa- 
tion]?] To take into account possible calibration issues 
between the datasets we added a multiplicative param- 
eter b which corresponds to the gain ratio and an ad- 
ditive parameter a, that compensates for possible DC 
offsets between the two datasets. We then fitted I V2 — 
b x f(T e ) + a with HT e ) given by Equation]! The bot- 
tom panel of Figurefin shows the observed 1 .5 GHz (I U2 ) 
vs 21.5 GHz (I Vl ) relationship together with the modeled 
function for the best fit T e = 1 1376 ± 1050 K. Note that 
in addition to the very low scatter about the observed re- 
lationship, the gain ratio (b) and DC offset (a) between 
the two datasets is 1.08 ±0.02 and 1.4± 1.0 MJy sr" 1 , re- 
spectively, which lie within the calibration uncertainties 
given in Table I] 

As shown in Figure [8] the relationship between free- free 
intensities involving 8.4 GHz (A-band) data are notice- 
ably scattered. This could be related to a flux calibration 
discrepancy between the GBT and VLA 8.4 GHz data. 
The overall relationship for 8.4 GHz vs 21.5 GHz inten- 
sities is linear as expected between two optically thin 
frequencies. Shown on this plot is the line with a fitted 
slope of 1.19±0.12. The predicted slope is 1.12 (assum- 
ing our T e = 11376K estimate). As the slope depends 
only weakly on T e , the difference between the measured 
and the expected values is plausibly related to the gain 
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Fig. 8. — Relationship between 21.5 GHz and 8.4 GHz (top panel) 
and 8.4 GHz and 1.5 GHz (bottom panel) free-free emission intensi- 
ties. For the bottom panel, we show the histogram of points density 
in gray-scale. Gray lines are the fitted relationship (see Section[7l. 

ratio between the two datasets and is consistent with cal- 
ibration uncertainties given in Table [l] The 1.5 GHz vs 
8.4 GHz intensity relationship (bottom panel of Figure [8| 
is non linear as predicted by Eq. [4] The best fit electron 
temperature in this case is T e = 13457 ± 3000 K, in good 
agreement with that from Figure]?] Due to the high scat- 
ter of the 1.5 GHz vs 8.4 GHz intensity relationship, we 
adopt the former (T e = 11376 ± 1050 K). 

This es timate is signifi c antly above the 7865 ± 360 K 
given by Subrahmanyan (19921 that was estimated di- 
rectly from the brightness temperature at the free-free 
emission peak in a 330 MHz map, where the emission 
should be optically thick. In M42 most of the ionized 
gas emission comes from the main ionizing front (MIF) 
located beyond 6' 1 OriC at the boundary with OMC1. 
This MIF is seen face-on. Therefore, when there is signifi- 
cant free-free optical depth there is a "photosphere" soon 
reached on the near side of the MIF. This is certainly the 



case for the measurement by Subrahmanyan ( 1992] ) . It is 
well established that the electron temperature rises from 
the ionizing star to the ionizing front, due to the harden- 
ing of the radiation field as the lowest energy photons are 
absorbed first. Usin g different ions emitt ing at succes- 
sively greater depths, |Esteban et al.| ( |1998[ ) has observed 
the increase of T e as one approaches the MIF on a single 
line of sight: T = 8730 ± 800 K from the Balmer lines- 
to-continuum ratio, T e [Ora] = 8300 ± 210 K, T e [Nn] = 
9850 ± 375K and T e [Sm]) = 103Q0±|^°K. Moreover, 
according to the modeling of the MIF with the software 
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Fig. 9. — Averaged spectrum of Orion KL/BN in t he 80" beam 
of ISO LWS. Triangles are ISO LWS data taken from |Lerate et al~| 
l |2006[ l. The dot-dashed line is the free- free emission component, 
tne clashed line is the thermal dust emission component and solid 
line is the total. For comparison, the dotted line shows the free- 
free emission without taking into account the opacity. This fitted 
function is described in the Appendix. 



packa ge Cloudy (08.00, latest described in Ferland et al. 
19981, the 330 MHz photosphere is located at the same 
depth into the Hn region as the [Om] emission. Since 
the opacity at 1.5 GHz (L-band) is t < 1 over most of 
the map, with our correlation method, we estimated the 
electron temperature weighted by the EM. The high- 
est densities and hence the highest opacities are at the 
neutral edge of the MIF so our estimate of T e is biased 
toward this region. 

8. DUST EMISSION 

Toward M42-OMC1, dust emission from (1) the cold 
dense part of OMC1, (2) the hot interface of the molecu- 
lar cloud heated by the Trapezium stars, and (3) the H n 
region, are mixed. In their study based on PRONAOS 
submillimeter o bservations and ancil lary data as short as 
200 and 90 /mi, |Dupac et al.| ( |2001[ ) derived a tempera- 
ture of 66 K for the dust. Hot dust fills a larger fraction 
of 3.5' PRONAOS beam than cold dust. Consequently, 
even if only from the H n region and the surface of OMC1 
(which has a relatively small column density compared 
to that of the cold component extending through the 
molecular cloud) , emission from hot dust is energetically 
important. This is the situation expected as the depth of 
the hot dust in OMC1 is set by when the opacity to UV 
light, tjjvi is approximately 1. Moreover, for the short- 
est wavelengths they used (90 an d 200 /mi), where the 



hot dust dominates the emission. Werner et al. ( 1976 1 
derived a dust temperature ranging between 55 and 85 K 
from the 50/100 /mi emission ratio in a 1' beam. 

At smaller angular scales toward the regions of the 
highest column density such as Orion KL, the colder dust 
component from the more shielded part of the molecular 
cloud must dominate. Especially at wavelengths longer 
than 350 /im, we expect to find a lower apparent dust 
temperature (even allowing for embedded energy sources 
such as BN and Ircl). Figure [9] displays the averaged 
spectrum of Orion KL/BN in the 80" beam of the In- 
frared Space Observatory long wavelength spectrometer 
(ISO LWS), using our dataset and three points at 80, 



et al. (2006). Fitting the function described in the ap- 
pendix, we obtained a dust temperature Td — 42 ± 3 K 
and a spectral emissivity index (3d = 1-3 ± .1 for the 
dust, lower than the average over the nebula. |Lis et al] 
(I1998P derived a value of (3 d = 1.8 toward Orion KL/BN 
from the 350/1100 /tin emission ratio and assuming a 
temperature of 55 K. They would have obtained an even 
higher value of (3d by using Td = 42 K. This difference 
could be related to their use of the 3300 /mi (90 GHz) 
data of |Salter et al. (19891 to subtract the free- free con- 
tribution from their 1100 /mi data. As we discussed in 
Section [6] and is evident in Figure [9] dust emits signifi- 
cantly at 3300 /mi (90 GHz) toward Orion KL/BN so that 
they must have overestimated the free-free contribution 
at 1100 /mi leading to an underestimation of the dust 
emission and then to an overestimation of (3d using the 
350/1100 /im emission ratio. 

Our dataset, summarized in Table [T] does not go out to 
short enough wavelengths to allow us to constrain both 
(3d and T d at other points. Therefore, to study (3d vari- 
ations between the free-free emission peak, Orion S and 
Orion KL/BN locations (shown on Figure [3| we fixed T d 
to the value of 42 K derived above. Spectra for these po- 
sitions are shown in Figure [5] and were obtained by aver- 
aging the signal within a 30 diameter circle. We derived 
values of (3 d ~ 2.0±0.1, 1.4±0.1 and 1.2±0.1 for the free- 
free emission peaks, Orion S and Orion KL/BN, respec- 
tively. At the Orion KL/BN location, (3d goes from 1.3 to 
1.2 as the area averaged over is decreased from 80" to 30" 
in diameter. The decrease of (3d toward den ser regions 
of OMC1 has been reported previously (e.g ., iGoldsmith 



et al.|1997|[Lis et al.|1998[|Johnstone fc Ba lly 1999). This 
may not be a physical effect, but rather the result of de- 
scribing the spectral energy distribution from dust that 
has a distribution of temperatures using a single value of 

In order to study the trend of Td related to the ob- 
served spectral evolution, we fitted the same three spec- 
tra with fixed (3 d = 1-8 (the low value of (3 d = 1.3 
obtained toward Orion KL/BN would give an unphysi- 
cally high temperature for the free-free emission peak). 
We got T d = 82 ± 10, 21 ± 1 and 15 ± 1 K for the free- 
free emission peaks, Orion S and Orion KL/BN, respec- 
tively. One must note that the obtained temperature 
at the free-free emission pe ak is then clos e to th e maxi- 
mum value T d ~ 85 K from |Werner et al.| ( [1976]) . More- 
over, Td ~ 82 K and (3d ~ 1.8 lead to peak emission at 
A = 30 — 40 jim. This is consi stent with the prediction 



100 and 150 /im taken from the ISO LWS data of|Lerate 



of the dust model described in Compiegne et al. ( 2008 ) 
with an intensity of the UV radiation field of x ~ 10 4 j 
compatible with the expected value at the hot interface 
of the molecular cloud heated by the Trapezium stars 
(x ^ 10 5 )- Toward OrionKL/BN, T d ~ 15 K is a lower 
limit ((3d ~ 1.3 from the spectrum of Figure [9] leads 
to Td ~ 30 K) but is nevertheless feasible for dust in a 
dense shielded environment. 

Finally, it is worth emphasizing that a spectral de- 
pendence on evolution of both temperature and spectral 
emissivity index of the dust is plausible. We have al- 
ready mentioned a decrease in temperature from ~ 85 K 
to ~ 20 K going from the highly excited environment 
around the Trapezium (the free-free emission peak) to 
the denser regions shielded from the UV radiation field 
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(Orion KL/BN and Orion S spectra). A decrease of the 
spectral emissivity index toward a denser medium can 
arise from dust grain properties evolution (shape, size, 
porosity, composition) following processes such as coag- 
ulation or accretion. It is therefore remarkable that al- 
though the observed spectra are averages along the lines 
of sight among different dust emission properties, a sin- 
gle dust temperature and spectral emissivity index suf- 
fice to reproduce the sp ectrum between 350 and 3300 /im. 
Goldsmith et al. ( 1997 1 also reported no systematic trend 
indicative ot frequency dependence of (3d between 450 /im 
and 1100 /im toward this region. 

9. CONCLUSION 

We have presented one of the first maps made with the 
MUSTANG camera that allows observations at 90 GHz 
using the 100 m diameter Green Bank Telescope (GBT). 
We obtained a 5' x 9' continuum map in the bright Huy- 
gens region of M42 at 9" resolution with a sensitivity 
of 2.8mJy/beam. Dual data analysis pipelines produced 
similar maps which are free from artifacts and have dy- 
namic ranges of over 1000. 

To analyze the physics of this region, we used a multi- 
frequency dataset ranging from 1.5 GHz to 854 GHz 
whose characteristics are summarized in Tabled] We 
also presented a newly obtained 150 GHz map from the 
GISMO camera on the 30 m IRAM telescope. 90 GHz is 
an interesting transition frequency, such that MUSTANG 
detects both the free-free emission characteristic of the 
H ii region created by the Trapezium stars, normally seen 
at lower frequencies, and thermal dust emission from the 
background OMC1 molecular cloud, normally mapped 
at higher frequencies. Using an analytical spectral func- 
tion (described in the appendix) including free-free and 
thermal dust emission, we studied these two components. 

We derived a peak value of the emission measure (EM) 
of 7.6 ± 0.8 x 10 6 pccm" 6 in a 33.5" FWHM beam (the 
resolution of our 21.5 GHz map). Using only higher res- 
olution radio frequency maps (~ 8.4"), we were able to 
get the spatial variations of the EM at the same resolu- 
tion as optical estimates by Baldwin et al. (19911. Both 
the absolute value and the spatial variation of our EM 
estimate are in close agreement with those of Baldwin 
etaLldlWl. 



We derived a new estimate of the EM averaged elec- 
tron temperature of T e = 11376 ± 1050 if. Our method 
is based on the fitting of the non-linear behavior of the 
relationship of free-free emission intensities at optically 
thin and optically thick frequencies whose only free pa- 
rameter is the electron temperature. Our temperature 
estimate is in agreement with previous estimates for the 
main ionization front T t ([S in]) 



103001*8° 



K (see 



Es- 



teban et al. 1998 1 where the density rises and is then 
expected to dominate the EM. 

We used ISO-LWS data in addition to our dataset in 
order to constrain both the temperature and spectral 
emissivity index of the dust within the 80" ISO-LWS 
beam toward Orion KL/BN. We derived T d = 42±3K 
and (3d = 1.3±0.1. By comparing spectra from the free- 
free emission peaks, Orion S and Orion KL/BN we have 
shown that Td and (3 d decrease from the Hn region and 
the excited OMC1 interface to the denser UV shielded 
part of OMC1 (OrionKL/BN, OrionS). While the tem- 
perature tendency is obviously related to the UV exciting 
radiation field and its extinction in the molecular cloud, 
the Pd tendency could be related to dust grain proper- 
ties evolution (shape, size, porosity, composition) follow- 
ing processes such as coagulation or accretion. Since the 
MUSTANG data presented in this paper was taken there 
have been significant improvements to MUSTANG'S sen- 
sitivity and more are expected in the future. MUSTANG 
is available to the community and we expect a bright fu- 
ture for 90 GHz science using the GBT. 
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APPENDIX 
FITTED SPECTRAL FUNCTION 



The function used to fit the spectra is 



with l Vi d the thermal dust emission intensity and I v ,ff the free- free emission intensity. We have 

h,d = C d B v (T d ,v) zA, 



(Al) 



(A2) 



with v the frequency, Cd a parameter relating to the dust opacity and column density, B v (T d , v) the Planck function 
at the temperature Td of the dust, and (3d the spectral index for the dust emissivity. The free-free emission intensity 
is expressed as 



Iujf = I Q<v ,ff{EM,T e ,v) x 



(1 



-t(EM,T b ,v) 



with 



t(EM, T e , v) 

I 0iV<ff (EM,T e ,is) = C ff T- 1 ' 2 g ff (T e ,u) e~& EM 



(A3) 
(A4) 
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being the intensity without opacity at frequency v for a value EM of the emission measure, an electron temperature 
T e and with C// = 5.4 10~ 39 taking the mean ion charge to be ~ 1 and the refractive index of the plasma ~ 1. We 
use the Gaunt factor approximation 



v/3 



ln(5 x 10 7 

TT V 



-■ill 



The optical depth is defined as 

t(EM, T e , v) 



.asxurVg)—^)- 



EM 
'pc cm - 



(; 



;)a(v,T e ) 



(A5) 



(A6) 



with a(u, T e ) the correction factor tabulated by Mezger & Henderson ( 1967[ ) that is close to 1. Among the 5 parameters 
{Cd, Td, (3d, EM, T e ) to fit, we fixed T e to 11376K, the mean value that we derive in Sectio n [7] (the dust parameters 
are not very sensitive to this choice). The fitting was done using MPFIT ( |Markwardt|[2009l ). 
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